1 Introduction

The goal of this project is to develop a robust and transparent framework for assessing afforestation suitability in Iceland. Afforestation plays a key role in climate adaptation, land restoration, and biodiversity enhancement, but its success depends on finding locations where environmental conditions, social priorities, and policy goals align. This project aims to identify such areas through an integrated, data-driven modern approach that combines ecological modelling, remote sensing, and spatial decision analysis.

The workflow begins with a broad-scale ecological assessment that identifies general patterns of forest suitability across Iceland. This first stage captures large-scale influences such as climate and soil, providing a foundation for understanding where afforestation is biophysically feasible. A second, finer-scale model then refines these predictions using high-resolution satellite imagery and terrain data, capturing local variations in topography, microclimate, and land surface characteristics. Together, these models provide a multi-scale understanding of where trees are most likely to establish and thrive.

Building on these ecological layers, the analysis incorporates social and infrastructural factors such as accessibility and proximity to population centres. These elements ensure that ecologically suitable areas are also practical for implementation, management, and community engagement. The inclusion of such human dimensions reflects the reality that successful afforestation must balance environmental potential with social and logistical feasibility.

Finally, the project applies a probabilistic Multi-Criteria Decision Analysis (MCDA)(suitability ananlysis tool in ESRI) to integrate all dimensions of suitability—ecological, social, and infrastructural—into a unified framework. Unlike traditional MCDA approaches that rely on fixed weights, this probabilistic method samples from plausible weight ranges to represent uncertainty in the weighting process. This produces a spectrum of possible outcomes rather than a single deterministic map, supporting more flexible, risk-aware planning.

Overall, this framework provides a transparent and scientifically grounded foundation for afforestation planning in Iceland. By linking large-scale ecological modelling with fine-scale environmental data and social considerations, it enables policymakers and land managers to identify where afforestation can deliver the greatest ecological and societal benefits.

2 Modelling workflow

3 Ecological suitability

The ecological suitability assessment forms the foundation of the afforestation suitability analysis workflow. We employ a two-step modelling approach: an initial broad-scale assessment provides a general overview of suitable regions across Iceland, followed by a refinement stage that incorporates high-resolution environmental data (satellite and terrain) to capture local variability. This hierarchical framework ensures that both large-scale patterns and fine-scale environmental nuances are considered, resulting in more accurate and ecologically informed afforestation planning.

Broad-scale model → Local-scale refinement

3.1 Broad-scale forest suitability mapping

Ecological niche modelling

To identify ecologically suitable areas for afforestation in Iceland, we first employ a broad-scale species ecological niche model (ENM). ENMs are a tool used to predict where conditions are environmentally suitable to support species survival and reproduction. By relating known occurrences of a species (such as observation or survey data) to environmental variables like temperature, precipitation, soil type, or elevation, ENMs identify the combination of conditions that define the species’ ecological niche. These relationships can then be projected across space to estimate potential distributions—such as identifying suitable habitats in unsurveyed areas. ENMs are widely used in ecology, conservation planning, and biodiversity management to understand species–environment relationships and guide decision-making (Shilky et al. 2023).

Species occurrence data

Here, we used the Natturulegt birkilendi and the raektad skoglendi datasets to gain an understanding about the current distributions of naturally occuring Birch and cultivated forests in Iceland. Since this data is given as polygons, 5000 points were randomly sampled from the polygons for subsequent ENM. In practice, rather than combining all species together as done here, species-specific ENMs could be used to obtain individual species preferential habitats across Iceland. Future work could focus on this.

Environmental Predictor Variables

The ENM incorporated climate, soil, and topographic variables derived from multiple global datasets. Climate variables were obtained from WorldClim, soil properties from SoilGrids, and terrain attributes from the ArcticDEM and are shown below.

Table: Environmental predictor variables used in the ecological niche model.
Dataset Description Units
WorldClim All 19 bioclimatic variables (temperature and precipitation predictors) —
SoilGrids Bulk density of the fine earth fraction cg/cm³
SoilGrids Cation Exchange Capacity of the soil mmol(c)/kg
SoilGrids Volumetric fraction of coarse fragments (>2 mm) cm³/dm³ (vol‰)
SoilGrids Proportion of clay particles (<0.002 mm) g/kg
SoilGrids Total nitrogen (N) cg/kg
SoilGrids Soil pH pHx10
SoilGrids Proportion of sand particles (>0.05 mm) g/kg
SoilGrids Proportion of silt particles (0.002–0.05 mm) g/kg
SoilGrids Soil organic carbon content in the fine earth fraction dg/kg
SoilGrids Organic carbon density hg/dm³
SoilGrids Organic carbon stocks t/ha
ArcticDEM Terrain slope degrees
ArcticDEM Topographic Position Index —

Broad-scale Ecological suitability Predictions

ENM predictions were estimated across Iceland at a 1 km spatial resolution. These predictions highlight areas with the highest ecological suitability for afforestation, providing a foundation for more detailed, fine-scale planning. The predictions can be seen below, where 0 indicates low suitability for tree establishment and 1 indicates highly preferential.

Figure 1. Broad-scale predictions of ecologically preferential areas for afforestation

3.2 Refining Local-scale suitability estimates

The next step is to make the broad 1km predictions more detailed. We do this using a method called residual stacking. First, we create a model at a much finer 10m resolution, using detailed elevation data and satellite images. This model picks up local differences that the coarse model misses. We then combine this fine-scale information with the broad-scale predictions to produce a more detailed map of areas suitable for forest growth.

Digital elevation models

We use the Arctic DEM to derive fine-scale topographic variables that influence forest establishment These include slope, elevation, topographic wetness index, and the amount of solar radiation received during the summer growing season. By incorporating these high-resolution elevation-derived features, the model can capture local variations in microclimate and soil moisture that the coarse 1km predictions may miss.

Satellite imagery

We used Google Earth satellite embeddings, which are computer-generated summaries of the Earth’s surface derived from multispectral satellite imagery. These embeddings are produced at a 10-metre spatial resolution and represented as a 64-layer stack, where each layer encodes a different aspect of the landscape. Rather than showing colours or images directly, these 64 layers form an abstract representation of what the satellite observes — including information both within and beyond the visible spectrum (such as infrared signals related to vegetation and surface materials). In simple terms, each embedding captures a compact “fingerprint” of the Earth’s surface, describing texture, reflectance, and land cover characteristics in numerical form. This allows us to analyse and compare locations consistently without working directly with the raw satellite imagery.

Table: Fine-scale environmental predictor variables used in the ecological niche model.
Dataset Description Units
ArcticDEM Terrain slope degrees
ArcticDEM Elevation (height) meters
ArcticDEM Topographic Wetness Index (TWI) unitless
ArcticDEM Mean radiation received during growing period MJ/m²
AlphaEarth 64-layered satellite embeddings unitless

Broad-scale Ecological suitability Predictions - Reykjanes peninsula

# Set layout to 2 rows, 1 column
par(mfrow = c(2, 1))

plot(broad_scale_cropped, main = "Broad-scale model")
plot(fine_scale, main = "local-scale refined model")

# Reset to default
par(mfrow = c(1, 1))

4 Social Factors

In our suitability analysis, social factors play a critical role in determining the optimal locations for afforestation. Two key social variables are considered: distance from roads and distance from towns. These factors account for the logistics of planting and maintenance, as well as public accessibility and recreational value, reflecting the preference for having forested areas within reach of nearby communities.

4.1 Distance from Roads

Accessibility to existing infrastructure is modeled using the inverse exponential function:

\[ P(x) = e^{-k \cdot x} \]

where \(x\) is the distance from the nearest road and \(k\) is a decay parameter controlling how quickly preference decreases with distance. Locations closer to roads are preferred because they facilitate management, transportation, and access for planting and maintenance activities.

This formulation reflects a diminishing preference with increasing distance — meaning accessibility rapidly declines as areas become more remote. A higher value of k would make the preference drop off more steeply, emphasizing proximity to infrastructure, while a smaller k creates a more gradual decline.

By expressing accessibility as a continuous exponential decay, this model ensures smooth transitions in suitability across space, avoiding abrupt thresholds that might otherwise create artificial boundaries in the spatial decision-making process.

4.2 Population-weighted Proximity Preference

We model afforestation preference as an exponentially decaying function of distance from towns, weighted by population size.
The preference at a location \(x\) is calculated as:

\[ P(x) = \sum_{i=1}^{n} w_i \, e^{-\lambda \, d_i(x)} \]

where

  • \(P(x)\) – preference at location \(x\)
  • \(w_i\) – population weight of town \(i\) (e.g., number of inhabitants)
  • \(\lambda\) – decay-rate parameter controlling how quickly preference decreases with distance
  • \(d_i(x)\) – Euclidean distance from location \(x\) to town \(i\)
  • \(n\) – total number of towns

This formulation captures how social accessibility to afforestation areas tends to decrease with increasing remoteness from population centers. Areas closer to towns — especially larger ones — are assigned higher suitability because they offer better community engagement opportunities.

The exponential decay ensures a smooth and continuous decline in preference with distance, avoiding abrupt changes or artificial thresholds. The weighting by population size (\(w_i\)) ensures that proximity to major towns contributes more strongly to overall suitability than proximity to smaller settlements.

In practice, the decay parameter (\(\lambda\)) can be tuned to reflect different planning priorities: a higher \(\lambda\) emphasizes planting near urban areas (steeper decline), while a smaller \(\lambda\) reflects a more gradual influence of distance, extending the social preference gradient further into rural landscapes.

5 Afforestation suitability mapping

Typically, in suitability analysis, a Multi-Criteria Decision Analysis (MCDA) framework is adopted to guide decision-making when there are multiple factors influencing the outcome. This approach underpins the Suitability Modeler tool in ESRI ArcGIS Pro, which combines various spatial datasets to identify the most suitable locations for a specific purpose. Users provide a set of weights that represent the relative importance of each influential factor — for example, environmental conditions, accessibility, or land characteristics. These weighted factors are then standardised and combined to produce a composite suitability surface, where higher values indicate areas that best meet the defined criteria.

In the current project, we extend upon traditional MCDA methods by introducing uncertainty into the weighting process. Rather than treating the weights for each factor as fixed values, we define them as plausible ranges that reflect uncertainty or variability in expert judgment and data quality. This approach recognises that the relative importance of each factor may not be precisely known and can vary under different assumptions or scenarios. By sampling across these ranges, we generate multiple realisations of the suitability model to explore how changes in the input weights influence the final results. This allows us to assess the robustness and sensitivity of the suitability outcomes, providing a more realistic representation of spatial uncertainty in decision-making. Probabilistic MCDA approaches are increasingly popular in environmental and spatial planning applications (Tang et al. 2018).

5.1 Normalization of criteria

Before performing the analysis, each raster layer is normalized to a common [0 to 1] scale to ensure comparability across different measurement units and value ranges. The ecological suitability model is already constrained between 0 and 1, whereas a min-max normalisation procedure is performed on the (a) Population-weighted Proximity Preference and the (b) Distance from Roads factors to re-scale the factors. This is a common re-scaling approach: \((x - min(x)) / (max(x) - min(x))\), where \(x\) is the variable of interest. This transformation scales the values of x to a 0–1 range, preserving the relative differences between observations while standardising the scale across different variables.

5.2 Suitability analysis

Once the social factors were min-max normalised, we stacked all layers and performed the probabilistic MCDA over the Reykjanes Peninsula using 100 iterations and the following weight ranges listed in the table below. Over each iteration, a value is randomly sampled from each range and combined into a set of weights that sum to one. These weights are then applied to the normalized raster layers in the suitability stack, generating a unique MCDA outcome for that iteration. By repeating this process across all iterations, we produce a stack of suitability maps that collectively capture the uncertainty in the weighting scheme.

Table: Weight ranges for probabilistic MCDA over the Reykjanes Peninsula.
Criterion Min_Weight Max_Weight
Eco 0.4 0.6
Road Proximity 0.2 0.3
Town Proximity 0.2 0.3

Below are preliminary suitability results for the Reykjanes peninsula, here a higher value indicates greater suitability, where 1 is highly suitable and 0 is unsuitable.

plot(suitability, main = "MCDA Mean suitability analysis results", col = viridis(100, option = "A"))

References

Shilky, BSPC Kishore, Gajendra Kumar, Purabi Saikia, and Amit Kumar. 2023. “Application of Species Distribution Modeling for Conservation and Restoration of Forest Ecosystems.” In Ecosystem and Species Habitat Modeling for Conservation and Restoration, 249–64. Springer.
Tang, Zhongqian, Shanzhen Yi, Chunhua Wang, and Yangfan Xiao. 2018. “Incorporating Probabilistic Approach into Local Multi-Criteria Decision Analysis for Flood Susceptibility Assessment.” Stochastic Environmental Research and Risk Assessment 32 (3): 701–14.